Algorithms for Molecular Biology

Conventional short read DNA alignment algorithms, such as widely used BWA-MEM2, devote a substantial fraction of their execution time to seed generation. We present RosaSeed, a short read alignment algorithm, based on a configurable seeding framework, that significantly accelerates alignment while maintaining accuracy comparable to that of BWA-MEM2. RosaSeed replaces the BWA-MEM2 seeding kernel and directly supplies candidate seeds to the existing BWA-MEM2 chaining and alignment extension pipeline, producing standard SAM output suitable for downstream variant analysis. RosaSeed was evaluated on both simulated and real short read datasets, and its seeding time, fastq-to-SAM time, alignment accuracy, and memory footprint were compared against BWA-MEM2, Minimap2, Bowtie2, and the Enumerated Radix Tree (ERT) algorithm. We also compared against ERT2, a proposed simplified variant of ERT that sacrifices a small amount of alignment accuracy to significantly increase the alignment throughput. With single-threaded software execution, our recommended RosaSeed configuration achieves 13.06x and 3.85x speedups in seed processing and fastq-to-SAM time, respectively, over BWA-MEM2. Corresponding seed processing (and fastq-to-SAM) speedups are 3.96x (2.48x) over ERT, 1.2x (1.44x) over ERT2, 4.82x (3.31x) over Minimap2, and 33.24x (10.83x) over Bowtie2. RosaSeed requires only 49.61 GB of peak memory, approximately 25% less than ERT and ERT2, while maintaining alignment accuracy comparable to that of BWA-MEM2. Following the recent release of minibwa, a prefetch-optimized FM index aligner, we additionally benchmarked RosaSeed against it on an upgraded 24-core, 48-thread AMD Zen 3 workstation. Our speed-optimized miniRosaSeed configuration achieves a 2.12x fastq-to-SAM speedup over minibwa for single-threaded execution while producing 2.58% higher standard accuracy. Against Strobealign, miniRosaSeed reported a 2.18x fastq-to-SAM speedup and 3.1% higher standard accuracy. For multi-threaded execution, miniRosaSeed was both faster and more accurate than minibwa for up to 26 software threads on the Zen 3 workstation. RosaSeed provides multiple runtime-memory configurations, spanning memory-efficient and high-performance strategies, enabling users to select operating points appropriate for diverse computational environments while maintaining alignment quality comparable to that of BWA-MEM2.

Multiple sequence alignments (MSAs) are often treated as a solved problem. However, ambiguously or incorrectly aligned regions are the source of much error in phylogenetics and comparative genomics. CRAIC is a cross-platform desktop workbench that treats an alignment as a set of hypotheses to be stress-tested. CRAIC uses the posterior match probabilities of a three-state pair hidden Markov model (pair-HMM) to derive complementary views of alignment confidence. The five views are (i) a reliability overlay with live column and residue masking; (ii) a multi-aligner disagreement map, (iii) a local realignment sandbox, (iv) an interactive posterior explorer, and, where the true alignment is known, (v) a per-residue view of correctly predicted homologies. When present, MAFFT, MUSCLE, ProbCons, Clustal Omega, ClustalW and PRANK are integrated in the tool. In benchmarks the reliability score can locate misaligned columns, but while CRAIC can implement automatic filtering, it is primarily built for inspection and curation.

We present a lossless compression algorithm for indexing a protein database while supporting fast taxonomic classification in the method Centrifuger. The algorithm is a new scheme of the previously proposed run-block compression algorithm to reduce the size of the Ferragina-Manzini (FM) index, and it scales better with alphabet size than the original. On the RefSeq prokaryotic and viral protein sequences, Centrifuger reduces the memory footprint by over a third compared to the method Kaiju that builds on a plain FM-index, while having comparable running time. Furthermore, the compressed FM-index is lossless and can locate matches of arbitrary length, which helps Centrifuger achieve greater accuracy than Kraken2, a k-mer-based taxonomic classification method. We leverage the computational efficiency of Centrifuger to create an index of size 182 GB for classifying the reads against the full nr database that contains about 250 billion amino acid characters. Using this index, Centrifuger reveals different SARS-CoV-2 infection states and viral transcriptome profiles across human cell types from single-cell RNA-seq data.

This study introduces a novel model that integrates the capacity of Potts models to capture structural dependencies with the computational efficiency afforded by Felsenstein's algorithm in phylogenetics. While the present work focuses on protein models, the PhyloPotts framework is equally applicable to alignments of discrete states, including RNA, DNA, and 3DI characters. Our findings indicate that, although phylogenetic estimates do not consistently exhibit reduced bias relative to standard models, the approach substantially enhances the accuracy of ancestral sequence reconstruction. The method is freely available via an open source implementation in the PhyloPotts package for BEAST 2 from https://github.com/rbouckaert/potts.

Nucleic acid secondary structure thermodynamic models have great utility for the analysis and design of nucleic acid structures, devices, and systems, yet their complexity and piecemeal development over decades pose challenges to further innovation. Employing different functional forms for each loop type within a structure, existing nearest neighbor (NN) models necessitate numerous loop-specific dynamic programming recursions in order to calculate physical quantities of interest. Here, we derive and implement a new class of models based on transfer matrices, yielding a single functional form for the free energy of all loop types (as the trace of a matrix product) and a single dynamic programming recursion for calculating physical quantities of interest. This new class of models has additional attractive properties: 1) incorporating a subensemble of states describing proximal interactions within each loop, 2) approximating logarithmic entropic scaling for all closed loops, 3) offering the opportunity to model correlations between the subensembles of neighboring loops by generalizing transfer matrices to tensors, and 4) providing a basis for systematically adjusting the number of model parameters to enhance accuracy as new empirical data become available for parameter regression. We use literature experimental data to regress parameters for representative new transfer matrix models, observing accuracy as good as or better than that of existing NN models, yet using dramatically fewer parameters (e.g., dozens vs hundreds or thousands). For complex ensembles containing an arbitrary number of interacting nucleic acid strands totaling N nucleotides, the new single recursion leads to exact dynamic programs for calculating diverse physical quantities (including the partition function, equilibrium base-pairing probabilities, minimum free energy and suboptimal proxy structures, Boltzmann-sampled secondary structures, and gradient of the partition function) with O(N3) time complexity. The simplicity of the new partition function dynamic program enables overflow-safe partition function calculations using renormalization to achieve scalability for large complexes (e.g., 50,000 nt) as well as formulation using a global matrix multiplication primitive that facilitates GPU acceleration, yielding speedups of an order of magnitude relative to NN partition function algorithms for large complexes. The new class of transfer matrix and tensor models provides an elegant framework for future efforts to advance the accuracy, speed, and scope of computational nucleic acid analysis and design.

For ten years, baltic (Backronymed Adaptable Lightweight Tree vIsualization Code) has been used to make annotated phylogeny figures in molecular epidemiology, including during the West African Ebola epidemic, the Zika epidemic in the Americas, and the SARS-CoV-2 pandemic, as well as other fields. At its core, baltic is a Python library used for the efficient parsing, traversal, manipulation, and visualisation of phylogenetic trees. It is a small library with few dependencies that reads tree formats common in phylodynamic analyses and gives the user leverage to interact with a lightweight tree data structure to produce publication-ready figures with matplotlib. We present baltic v1.0, its first formally released and documented version. baltic reads and writes BEAST Nexus, Newick, and Nextstrain/Auspice JSON, and can process large BEAST posterior tree files in parallel to extract user-defined posterior statistics. The same objects are used for tree manipulation and for plotting in a single script. New to this release: a set of rooting methods (midpoint rooting, rerooting on any branch, and root-to-tip regression); support for reticulate evolution, with reassortment and recombination edges; and composite figures that combine a tree with other data, such as Muller plots, skygrid plots, tanglegrams, and plots connecting trees to maps. The release includes a documentation site with an API reference, tutorials, and a matplotlib-style gallery of worked examples.

Summary: Ensembl VEP is the de facto reference for variant consequence annotation, but its Perl implementation limits throughput, and the open reimplementation benchmarked here is validated on a set too small to bound its disagreement with it. We present vep-rs, a Rust reimplementation of Ensembl VEP release 115.2, audited over 260,638,766 VEP consequence tuples, each a (location, allele, feature) key carrying a Sequence Ontology term set, from six variant datasets (ClinVar, gnomAD chr21, and 1000 Genomes chr21 across GRCh37 and GRCh38). On SNPs and indels vep-rs emits exactly VEP's tuple count on all six and attains F1 [≥] 0.999974, differing from VEP on 3,169 tuples, at most one in 38,000 per dataset; with two documented VEP defect shapes set aside, F1 rounds to 1.000000 on all six. Per consequence class, 29 of the 30 Sequence Ontology terms carrying at least 1,000 tuples exceed F1 0.99, and the lowest, 0.540 on start_retained_variant, is entirely a VEP defect. On structural variants, F1 is 0.975395 (GRCh37) and 0.909998 (GRCh38). Timed on 20 independent machines per engine and architecture, vep-rs is faster than VEP by 101x-284x on ARM (geometric mean 176x) and 78.3x-192x on x86 (geometric mean 135x). Five divergence classes are documented, four of them defects in VEP release 115 that vep-rs does not reproduce: a consequence call that contradicts itself, a splice-region term dropped from a variant that lies in a splice region, annotation against a chromosome the variant does not lie on, and output that depends on input batch composition. A pipeline keyed on VEP consequence terms can therefore adopt vep-rs against a measured, per-term bound, shed VEP's contradictory and batch-dependent calls, and annotate population-scale inputs at 78.3 to 284 times VEP's throughput, and per-file structural-variant sets at 22.5 to 52.9 times. Availability and implementation: vep-rs, its harness, both comparators, and every per-machine wall time are released under Apache-2.0 at https://github.com/natera-open-source/vep-rs (release v0.1.0, archived at https://doi.org/10.5281/zenodo.22837897). Supplementary information: Supplementary data are available at Bioinformatics online.

Recent work applies large language models (LLMs) to single-cell data, but the cell usually reaches the model as a short ranked list of genes. This list drops much of what defines a cell, so the model reasons from a partial view. Yet the information is not lost in the measurement, only in the text. Here, we introduce Cell-Lens, a training-free structured representation of the measured cell. It writes the cell as typed blocks, from paired measurements to marker-supported programs. On CellVerse, SOAR, CellPuzzles, and SC-Arena, across five model families, Cell-Lens improves performance by up to 27.1 points. It also narrows the model gap and reduces API-reported reasoning-token usage. The largest gains occur when a defining RNA marker falls out of the ranked list. In these cases, the paired protein data still capture the corresponding signal. Shuffling the added blocks removes the gain. The gain persists after selected label-associated proteins are removed. These controls link the improvement to cell-matched biological content. We also build and release CellSpectrum, a broader benchmark spanning eight sources and seven tasks with paired modalities. The gains hold there across tasks, tissues, and modalities. These results show that a better cell representation can improve single-cell reasoning without retraining the model. Code and CellSpectrum are publicly available at https://github.com/LiZaiyuan0619/Cell-Lens.

Motif scaffolding builds a protein around a prescribed structural motif -- an active site, a binding epitope, a metal site -- so that the motif is presented in its required geometry. The task is usually posed with the motif pinned down: its position along the sequence is given, and several rigid fragments inherit their mutual pose from a native structure. Many applications supply neither. We present NomadFlow, unindexed multi-motif scaffolding on all-atom SE(3) flow matching: the motif enters as a condition, a group projection preserves the internal geometry of each group exactly, and the groups move freely relative to one another until a set time, after which they act only as a condition. We assemble a preliminary benchmark for this setting, varying the granularity of the partition and the structural context around the motif, and obtain unique successes on 19 of 21 tasks. The trained model proposes a diversity of valid mutual poses, and on several tasks the native arrangement is among them. The shape of that solution space, however, varies strongly with the geometry of the condition: across tasks the successful layouts form a wide cloud, collapse onto the native pose, or contract to a connected low-dimensional family -- a continuum of poses rather than a set of isolated solutions. Multi-motif methods are therefore better compared by the space of layouts they propose and the successful part of it they retain than by yield or fold diversity alone. Code and weights: https://github.com/ipermyakov/nomadflow.

Motivation: Tree sets appear naturally when performing model based phylogenetic analyses. When tracking host-parasite co-evolution, testing sensitivity of different model assumptions, or even just checking convergence of individual runs of the same analysis, the question arises: are a given pair of tree sets the same? And if they are not the same, how do they differ? Results: We provide DensiTree 3, an interactive tool for inspecting pairs of sets of trees that allow visualising both tree sets using a single panel. Incorporation of CCDs provide better summary trees, taxon orderings and extra functionality like rogue detection, making it easier to interpret tree samples. Furthermore, a clade comparison graph allows inspection of differences in clade support and clade heights linked to the double-DensiTree panel. Together, these visualisations of tree sets make it easy to identify potential differences in the tree sets and interactively inspect these differences. Furthermore, recent advances in conditional clade distributions (CCDs) to represent tree distributions have been incorporated in DensiTree 3, and the most visible part of the user interface has been updated for a more modern look. Availability: DensiTree 3 is freely available from https://cs.auckland.ac.nz/~remco/DensiTree/. The program is licensed under GPL and source code is available via the link above.

Phylogenetic trees inferred from different genomic regions, or from organellar and nuclear genomes, frequently disagree owing to incomplete lineage sorting, hybridization, introgression or horizontal gene transfer. Such discordance is informative, yet the standard tool for comparing trees visually, the tanglegram, is restricted to two trees and is known to misrepresent their similarity. Here we introduce phylogenetic parallelograms, drawings in which any number of rooted trees are embedded into a common scaffold, a rooted phylogenetic network that displays all of them, and are drawn in parallel so that shared branches align while conflicting branches visibly diverge. On synthetic tree pairs, the complexity of a parallelogram tracks the topological distance between the trees, whereas the complexity of an optimized tanglegram does not. We apply the approach to genome-wide introgression in the Anopheles gambiae species complex and to organellar discordance in cats, grasses and Fagaceae. PhyloParallelograms is an open-source interactive application.

Integrative analysis of the properties of multiple networks pertaining to a biological system is a key problem in systems biology. But computing the properties of thousands of large networks poses a challenge. Most techniques that partially address this challenge, such as parallelism and optimized algorithms, treat the analysis of each network as an independent task. However, the networks arising from the same biological system are often similar to each other, and can we harness this similarity to achieve speedup? Towards this end, we propose a phylogeny-like data structure, nc-tree, to compactly represent an ensemble of related biological networks; and a versatile framework, ProcessNets, that runs incremental algorithms over the nc-tree to speedup property computations on the input networks. Our asymptotic analysis of space and time complexity, coupled with empirical evaluations of diverse ensembles of simulated and real-world GTEx gene coexpression networks, demonstrate that ProcessNets leads to significant space and time advantages when computing different network science measures on ensembles of similar networks. For instance, we achieve a compression factor of 2.1x and 8.1x in storing two representative ensembles of GTEx Muscle Skeletal networks with 1000 nodes/genes, and a time speedup of 2.86x to 3.97x over the fastest baseline for computing degree centrality in ensembles of networks with at least 14,500 nodes derived from subsampled GTEx datasets of ten tissues. These results are promising and encourage application of ProcessNets to analyze biological network ensembles derived from rapidly accumulating consortium/biobank-based datasets.

Motivation: De novo genome assembly relies on accurate path selection through complex graphs, yet unitig graph branching points remain a major ambiguity source. Multiple extensions may be equally plausible from topology alone, but sequence-level consistency of competing paths remains largely unexplored. Model: We introduce RASTRO, a reference-free framework for scoring candidate paths at unitig branching points using nucleotide sequence information. RASTRO-S employs a sample-specific n-gram model to estimate extension plausibility via bits-per-base surprisal, while RASTRO-L uses the Evo2 pretrained DNA language model for autoregressive path scoring. Both operate directly on unitig graphs without reference genomes or assembled contigs. To evaluate branch decisions independently, we develop a read-based validation strategy using boundary k-mers spanning the junction between graph context and candidate extensions. Results: Experiments on HIV-1 and SARS-CoV-2 unitig graphs from Logan show sequence-based models provide meaningful match with read-supported choices. Among 20,000 sampled branching points, unique read-supported choices were identified for 15,319 HIV-1 branches (76.6%) and 19,882 SARS-CoV-2 branches (99.4%). RASTRO-S improved match from 54.1% to 75.8% for HIV-1 and from 58.6% to 76.4% for SARS-CoV-2 as context length rose from n=5 to n=10. RASTRO-L also improved with context length, reaching 49.96% and 49.25% match. Both outperformed random selection and the longest-extension greedy baseline. Conclusion: RASTRO provides an efficient, interpretable framework for evaluating ambiguous path decisions in unitig graphs using raw-read evidence rather than contig-level alignment. Results demonstrate the potential of lightweight, sample-adaptive sequence models for improving local decision-making in genome assembly workflows.

Single-cell sequencing provides a rich source of data that can be used to model and estimate the evolution of somatic cells within individuals. Many common models for modeling DNA evolution in phylogenetics are based on continuous-time Markov processes that operate on aligned sites. For haploid DNA, such as found in mitochondria, there are many well established models for these processes but since nuclear DNA is diploid, containing copies of both the maternal and paternal genome, it is natural to develop continuous-time Markov processes adapted to diploid DNA. One possible model is that both the maternal and paternal evolve independently according to haploid models so that, in principle, DNA comes with a phase. i.e. a label that tells if the site was derived from the maternal or paternal sequence. However, this information may not be available in practice and only an unphased diploid genotype is available. That is, an unphased model groups both of the possible maternal-paternal labels into a single state. This grouping process, called ``lumping'', is used by other methods as well and results in a coarser process. Care must be taken so that the coarser process can still be modeled as a continuous-time Markov chain; the lumped process of the original phased Markov chain may or may not be Markovian, depending on restrictions imposed on the generator matrix. Here, we explore the mathematical formulations of phased and unphased Markov processes, as well as a collection of other diploid DNA models from the literature.

Selection-intensity estimates from RELAX are reported as a point value of K with a likelihood-ratio P. We report that, at default settings and on data of ordinary size, the program does not reproduce its own fits. Of 27 enzyme entries refitted under two optimiser configurations, none reproduced its log-likelihood to within 0.01 units; the median change was 103 units, the largest over 3,400, and four verdicts reversed. Eighty null orthologues reproduced none. A byte-identical command returned a distinct likelihood on every repetition, single-threaded, across three releases, and on alignments simulated under the fitted model, where 3.3 per cent of replicates reproduced. The documented random-number seed never reaches the generator when assigned on the command line, yet reads back as the value supplied. PAML localises the cause: its two-ratio model, without site classes, reproduced its log-likelihood for all 288 genes; its site-class models agreed for 27 to 67 per cent. The instability follows the mixture over sites, not the program. The output does not show it: 46 of 410 fits ended with a negative likelihood-ratio statistic, impossible under convergence, and 123 of 410 report a K re-estimated under a domain restriction rather than the unconstrained maximum. Of 234 published studies using RELAX, none reported a seed. Seeding while holding the thread count at one reproduced sixty of sixty runs on twenty genes under two releases; the seed alone reproduced none of five, and no documentation states the second condition. We recommend that fits be repeated and their dispersion published.

Genotype calling of sequencing data has improved over recent years thanks to technological and statistical breakthroughs. One such advance is the use of imputation-based genotyping methods in human genetics. For haploid and mostly non-recombining organisms, a natural equivalent is to leverage phylogenetic signal for imputation-based genotype calling. SARS-CoV-2 is one such case, and the use of improved genotyping methods in this species is particularly pertinent given the well-documented systematic problems with read alignment and assembly in rapidly evolving regions of its genome, motivating the need for statistically principled methods for genotype correction that complement existing community efforts to curate sequencing data. We present a statistical framework to combine phylogenetic prior probabilities with the information from aligned reads, and demonstrate its utility in correcting consensus calls made on short-read, viral sequencing data. We derive phylogenetically-informed posterior probabilities based on genome placement in a reference phylogeny, and read likelihoods calculated from machine-reported base quality scores. We further modify the prior at each site by estimating a maximum-likelihood scaling parameter jointly with reads at all other sites. We find that the phylogenetically-informed base calling method stays highly accurate across simulated depth and error conditions. In real datasets, the method vastly improves call rate while still maintaining similar accuracy as reference-based assembly methods. In problematic genomic regions, the phylogenetically-informed method confers both accuracy and call rate gains. Accuracy and call rate do not stratify with sequencing center and primer scheme used, suggesting broad applicability across sequencing technologies and potentially diverging sequences. Our results in real and simulated data demonstrate that phylogenetically-informed genotyping correction results in more complete sequences, achieving a higher call rate without compromising accuracy. In cases of low mean depth, sparse coverage, or difficult-to-sequence regions that may still retain a low call rate, the method improves the accuracy of genotype calls.

Shotgun sequencing of DNA mixtures has the potential to become an important tool in forensic genetics, particularly for degraded samples or complex mixtures. However, the analysis of shotgun-sequenced mixtures is challenged by several confounding factors, including sequencing errors, population structure, the potential presence of relatives, linkage, and linkage disequilibrium. Additionally, the statistical methods currently available for such data are limited. Here, we develop likelihood ratios that summarize the evidence for the presence of a suspects DNA in a mixture, assuming genotypes of the suspect and the victim are available, while accounting for DNA from unknown contributors modeled either as general population contamination or as a single contaminating individual. We show that these likelihood ratios are not robust to model assumptions and can falsely support the presence of the suspects DNA when the true contributor is even a distant relative. We therefore propose a combined evidence score, {Lambda}: the likelihood ratio minimized over contamination models, reported only when the evidence favors the suspect over a hypothesized sibling by at least a factor of C (default C = 10) under at least one contamination model, and otherwise set to 1 (log {Lambda} = 0). In simulations based on 1000 Genomes data, {Lambda} retains nearly all power to detect true contributors at mixture proportions as low as 5% while eliminating the false positives produced by model mis-specification, population structure, and by untested relatives. The methods are implemented in the open-source program DNAMIXTURE, providing the forensic and population genetics community a resource for robustly identifying small individual contributions to a DNA mixture.

Learning from imbalanced data remains a major challenge for graph neural networks (GNNs), as minority nodes are not only rare but also structurally marginalized within the graph. We address this issue with CLARA, a hierarchical learning framework that decomposes node classification into two stages: a coarse subgraph-level classifier that selects regions likely to contain minority instances, followed by a fine-grained node-level predictor within these regions. This design improves sensitivity while maintaining scalability, filtering out irrelevant areas and focusing learning on topologically meaningful neighborhoods. Experiments on benchmark graph datasets demonstrate substantial gains over established imbalance-handling methods, with CLARA reaching an F1-score of 88.3\%. The same strategy achieves significant improvements in protein-ligand binding site prediction, underscoring its broad and consistent effectiveness across both biological and general graph learning tasks.

Motivation: High-throughput sequencing technologies have driven rapid growth of biological sequence databases. Public repositories must therefore rely on automated computational heuristics to screen submitted sequences for errors and low quality. For example, SILVA SSU Ref, which exploits the conserved nature of 16S and 18S rRNA sequences, applies strict algorithmic quality controls yet still accepts sequences with up to 30% of their nucleotides deviating from any previously accepted sequence. This permissiveness creates opportunities for the admission of modified sequences, such as biologically plausible sequences generated by DNA foundation models. The vulnerability of public sequence databases to becoming polluted or poisoned with modified sequences necessitates the development of methods to detect such sequences. Results: We present the first investigation, to our knowledge, of the detectability of modified sequences. We consider simple computationally-modified 16S rRNA sequences that pass the quality control inclusion criteria of the SILVA database, which we generate via random substitutions. We present classifiers that can distinguish such modified sequences from natural 16S rRNA using conserved motifs. Our best classifier achieves over 90% sensitivity and specificity on our testing set when using a 5% artificial mutation rate. One feature used in our classifiers, gapped k-mers constructed from universally conserved nucleotides, was conserved across all three domains of life despite relying on exact matches to patterns found in E. coli, advancing our understanding of conserved grammatical structure in small subunit rRNA sequences. Availability and Implementation: Our source code is available at https://github.com/rainhaworth/16S-Mutation-Classifiers.

The genetic code is the rule that (apart from minor modifications) is used by all contemporary life on Earth to translate genes into proteins. It is summarized by a 4x4x4 table which maps the 64 possible triplet RNA codons to the 20 canonical amino acids ("AAs"). Understanding the origin of this code is a fundamental open problem in biology. We review the idea of a comma-free code, and put forward a general reason to suspect the primordial code was comma-free (even though the modern code is not). We explain why this primordial code should then persist as a comma-free subcode within the modern code, and identify the five optimal candidates for such a subcode. We then introduce the idea of "compatible comma-free codes" and explain how they allow multiple codes to evolve in parallel, circumventing a standard argument suggesting that code evolution should effectively freeze when the associated genome becomes sufficiently complex. In particular, we highlight the possibility that the original code was a comma-free doublet code, which was then overtaken by a compatible comma-free triplet code, which finally evolved (after the evolution of the start codon and ribosome) to the modern saturated triplet code. To assess this idea, we classify all inequivalent compatible comma-free doublet-triplet codes and, among these, identify the three most promising predecessors of the modern code. Consistency between these three best codes and the "wobble rules" in gene translation provide non-trivial evidence for this possibility.